Dynamical Instabilities in High-Obliquity Systems 

D. Tamayo 

Department of Astronomy, Cornell University, Ithaca, NY 14853 

O 

CN ■ dtamayo@astro.cornell.edu 
> 

o 

:z;: 

o 

J. A. Burns 



o 

in: 



(N 



and 



Ph " Department of Astronomy & College of Engineering, Cornell University, Ithaca, NY 14853 



and 

D. P. Hamilton 

Department of Astronomy, University of Maryland, College Park, MD 20742 



and 



> 

(N 
O 

O ! P. D. Nicholson 



Department of Astronomy, Cornell University, Ithaca, NY 14853 
Received ; accepted 



- 2 - 
ABSTRACT 

High-inclination circumplanetary orbits that are gravitationally perturbed by 
the central star can undergo Kozai oscillations — large-amplitude, coupled vari- 
ations in the orbital eccentricity and inclination. We first study how this ef- 
fect is modified by incorporating perturbations from the planetary oblateness. 
Tremaine et al. (2009) found that, for planets with obliquities > 68.875°, orbits 
in the equilibrium local Laplace plane are unstable to eccentricity perturbations 
over a finite radial range, and execute large-amplitude chaotic oscillations in ec- 
centricity and inclination. In the hope of making that treatment more easily 
understandable, we analyze the problem using orbital elements, confirming this 
threshold obliquity. Furthermore, we find that orbits inclined to the Laplace 
plane will be unstable over a broader radial range, and that such orbits can go 
unstable for obliquities less than 68.875°. Finally, we analyze the added effects 
of radiation pressure, which are important for dust grains and provide a natural 
mechanism for particle semimajor axes to sweep via Poynting-Robertson drag 
through any unstable range. We find that generally the effect persists; however, 
the unstable radial range is shifted and small retrograde particles can avoid the 
instability altogether. We argue that this is occurs because radiation pressure 
modifies the equilibrium Laplace plane. 

Subject headings: celestial mechanics — planets and satellites: dynamical evolution and 
stability 
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1. INTRODUCTION 

Satellites in inclined circumplanetary orbits that are subject to gravitational 
perturbations from the Sun can undergo large-amplitude eccentricity oscillations through 
the Kozai mechanism (Kozai 1962; Lidov 1962), often with catastrophic results (e.g., 
Carruba et al. 2002). As discussed below, such dramatic increases of eccentricity can occur 
as the pericenter slows its precession, allowing the solar tugs to systematically remove 
angular momentum from the orbit over part of the precession cycle. When the dominant 
perturbation is the Sun's gravity, this halting of the pericenter can only be achieved for 
highly inclined orbits (> 40°). In this paper we consider situations where additional 
perturbations are also important, thereby providing new ways to slow pericenter precession 
and to consequently generate large eccentricities. 

An important additional potential to consider is that due to the central planet's 
oblateness. The study of this perturbation's effect on satellites in combination with the 
Sun's gravity dates back to investigations of Saturn's moon lapetus by Laplace (1805), and 
later, by Tisserand (1896). Allan & Cook (1964) subsequently generalized this analysis to 
an arbitrary number of perturbers. These works were limited to circular satellite orbits, for 
which the motion can be expressed in terms of elementary functions. In the general case of 
eccentric orbits, however, the evolution is no longer integrable. In the special circumstance 
where the obliquity is zero, Kozai (1963), found a class of solutions where the argument of 
pericenter librates around ±90°, qualitatively similar to Kozai cycles. Lidov & Yarskaya 
(1974) tabulate and explore the integrable cases in the problem. 

Kudielka (1994) and Vashkov'yak (1996) built on these works and discovered solutions 
where most or all of the orbital elements remained stationary. However, they limited their 
analysis to low obliquities, applicable to the Earth-Moon system. Tremaine et al. (2009), 
henceforth TTN, analyzed the full range of obliquities and found the stationary solutions 
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for both circular and eccentric orbits. They further provide maps of the stabihty of these 
equihbria to eccentricity and angular momentum perturbations. Most importantly for 
this paper, they discovered that orbits around planets with obliquities > 68.875° undergo 
chaotic, large-amplitude oscillations in eccentricity and in inclination over the radial range 
from the planet where the two perturbations are comparable (see Fig. 9 in Tremaine et al. 
2009). For a visualization of the effect, see the orbital histories shown in the figures below. 

In this paper, we first investigate this case (including oblateness and gravitational solar 
perturbations) in a manner complementary to TTN. We use orbital elements in preference 
to TTN's vector approach, and we derive our results from the simple condition that the 
peri center be able to halt, rather than from the stability of the Laplace surface. We thereby 
sacrifice some mathematical elegance and generality in order to provide a more physically 
intuitive picture. In Section 3 we extend the work of TTN, which only considers orbits in 
the equilibrium plane, to consider the general case of orbits inclined to this Laplace plane. 

Although each added perturbation greatly increases the system's complexity, in Section 
4 we incorporate radiation pressure so as to be able to study the motion of dust grains, for 
which such forces matter (Burns et al. 2001). Most importantly, radiation forces generate 
Poynting-Robertson (P-R) drag. P-R drag causes an orbit's semimajor axis to decay 
(Burns et al. 1979), allowing it to sweep through the radial range from the planet in which 
the eccentricity becomes unstable. We investigate whether radiation pressure's additional 
effects alter this unstable radial range or are even capable of stabilizing particle orbits 
against the instability found by TTN. 
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2. EVOLUTION UNDER PERTURBATIONS FROM SOLAR GRAVITY 

AND PLANETARY OBLATENESS 

2.1. Kozai Oscillations 

We start by reviewing tlie features of tlie Kozai meclianism tliat are essential to our 
work in order to motivate the strategy pursued in the rest of the paper. In our context, 
Kozai oscillations result solely from the Sun's gravitational perturbations on a body in an 
inclined circumplanetary orbit. 

From a planetocentric perspective, the Sun "orbits" the planet in the latter's orbital 
plane; to avoid confusion with the particle's orbital plane, we will hereafter refer to the 
planet's orbital plane as the "ecliptic" (even though the latter term strictly refers to the 
Earth's orbital plane). When interested in secular timescales much longer than the planet's 
and particle's orbital periods, one can time-average over the Sun's and particle's orbits 
and treat their masses as distributions smeared over their paths in the sky. Furthermore, 
since a circumplanetary particle lies much closer to the planet than to the Sun, it is usually 
sufficient to take only the leading quadrupole term in an expansion of the solar potential 
in powers of a/a^, where a is the circumplanetary particle's semimajor axis, and is the 
planet's semimajor axis. We point out, however, that Katz et al. (2011), Naoz et al. (2011) 
and Lithwick & Naoz (2011) have found that including the octupole term can introduce 
qualitatively different phenomena including flips from prograde to retrograde orbits, and 
eccentricities arbitrarily close to unity. In the limit where the circumplanetary particle's 
mass is negligible, they find that the octupole correction can be ignored when eM 1, 
where 



and Cp is the planet's orbital eccentricity. In this paper we consider only the Sun's 
quadrupole potential, and our results are therefore only applicable to cases where e^/ ^ 1. 
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The secular problem truncated at quadrupole order was first analyzed by Kozai (1962) 
and Lidov (1962). We choose to work with the orbital elements (a, e, i, Q, u). The equations 
of motion in these variables are given by (Innanen et al. 1997; Carruba et al. 2002, though 
see an erratum common to both papers in Carruba et al. 2003): 



de 15eQn 2\i/2 • 2- • o /n\ 

— = — - — e(l - e ) ' simE sm2uE, (2) 
dt 8 



diE 15e0n 
It ~ 16~ 

dt 4 ^ ^ 



e^{l-e^)-^^^sm2uE sm2iE, (3) 



2(1 - e^) + 5sin^W£;(e^ - sin^^i?) 



(4) 



where n is the particle's mean motion, e is the particle orbit's eccentricity, ue its argument 
of pericenter, and tE its inclination to the ecliptic, in which the Sun moves. The subscript 
E has been added to the angular quantities to emphasize that they are measured relative 
to the ecliptic plane. The quantity characterizes the strength of the solar perturbation 
relative to the dominant planetary gravity and depends on the distance from the planet; it 
is given by 

'® M,a3(l-e2)3/2' 
where Mp and Mq are the planet's and Sun's masses, respectively. 

The fact that the Sun has been averaged over its orbit and that the potential it creates 
is therefore time-independent means that energy (and thus a) is conserved. Furthermore, 
Kozai (1962) realized that the problem's symmetry guaranteed the conservation of the 
component of angular momentum perpendicular to the ecliptic, Lz = -y/GMpa(l — e^) cosi^;, 
where G is the gravitational constant. This renders the system a one-degree-of-freedom, 
integrable system in the (e, ue) plane, i.e., one can divide Eq. 2 by Eq. 4, eliminate ie 
using Lz, and solve for e as a function of uje- For initial values of 6 = Vl — cos < 3/5 
{iE > 39.2° for e ^ 1), a stable equilibrium solution exists where e,iE, and ue are 



- 7- 



stationary. In this case, the phase portrait in the (e, cu^) plane consists of two types of 
solutions: 1) ones that trace out paths around the stationary point so that Ue librates 
between minimum and maximum values, and 2) ones where ue circulates. For O > 3/5, no 
stationary point exists, and only circulating solutions are possible. These behaviors can be 
seen in Figs. 2-8 of Kozai (1962) and Fig. 2 of Carruba et al. (2002). 

More qualitatively, Eq. 2 indicates that the pericenter's orientation within the orbital 
plane (given by ue) determines whether the eccentricity grows or shrinks. Normally ue 
circulates swiftly, as the term in brackets in Eq. 4 is roughly constant for small e and i. 
This results in a small-amplitude eccentricity oscillation (due to the sm2uE term in Eq. 
2). However, if duE/dt ever approaches zero in an orientation where sin2a;£; > 0, the 
eccentricity can grow to large values. This can occur in the Kozai problem whenever the 
relative inclination, iE, between the particle's orbit and the distant perturber's orbit is 
significant. For small eccentricities, Eq. 4 indicates that duE/dt equals zero for some values 
of ue when sin^iE > 2/5 (i.e., ie > 39.2°). 

Large inclinations to the ecliptic therefore provide one way for the eccentricity of 
circumplanetary orbits to grow to large values; however, adding other perturbations may 
allow for additional possibilities. 

2.2. Adding Planetary Oblateness (J2) 

A planet's oblateness, represented by its J2 coefficient, causes pericenter precession but 
does not produce secular effects on an orbit's eccentricity or inclination (e.g., Danby 1962). 
One can therefore imagine that orbital configurations may exist in which the ue precession 
from J2 cancels that from the Sun, making ue constant and allowing the eccentricity to 
grow to large values according to Eq. 2. 
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The conservation of mentioned in the previous section is due to the quadrupole 
potential's azimuthal symmetry, and this causes the eccentricity and inchnation evolution 
to be coupled (Kozai 1962). But when one adds planetary oblateness, which is invariant 
about a different axis (the planet's spin pole), this symmetry of the classical Kozai case 
is destroyed. Hence, the eccentricity and inclination become decoupled, and the optimal 
choice of a reference plane from which to measure all angles is no longer obvious. 

An appropriate choice is the local Laplace plane, which lies between the planet's 
equatorial plane and the ecliptic. If a particle on a circular orbit has its orbital plane 
align with this equilibrium plane, the torques from the Sun and J2 balance so that the 
orbit's angular momentum vector remains fixed, and the orbital plane does not precess 
(Allan & Cook 1964). A circular orbit not aligned with the local Laplace plane will have its 
orbital axis precess around the equilibrium Laplace plane axis (see Fig. 1). This represents 
a compromise between the particle orbit attempting to precess around both the planet's 
spin axis and the Sun's orbital axis. 

More generally, the orbit normals of eccentric orbits will wobble as they undergo their 
precession cycle since their changing distance from the planet means they "sense" a range 
of Laplace planes. We note that even in cases where the orbital plane itself does not precess, 
the pericenter may still precess within that orbital plane. It is specifically the ability of the 
pericenter to halt that gives rise to large eccentricities, as argued in the first paragraph of 
this section. 

Because the strengths of the two relevant perturbations vary differently with distance 
from the planet, the local Laplace plane shifts as the particle's semimajor axis varies. Near 
the planet, where the torques on the orbit are predominantly caused by oblateness, the 
Laplace plane nearly coincides with the planet's equatorial plane. Far from the planet, 
where solar torques dominate, the Laplace plane aligns closely with the ecliptic (in which the 



Fig. 1. — Normal to the local Laplace plane z lies between, and is coplanar with, the planet's 
spin pole fip and the ecliptic normal Uq. The normal to an arbitrary particle's orbit plane 
j will precess around z at approximately constant inclination i, sweeping out a cone. The 
obliquity (pQ is simply the angle between the vectors rip and n©, and (p represents the angle 
between rip and the z axis. As the semimajor axis changes and the relative strengths of the 
Sun's and planet's perturbations vary, the Laplace plane will shift, and will vary. 

Sun "moves"). Between these limits, the Laplace plane takes on intermediate orientations, 
generating a warped Laplace surface. The Laplace plane at a given semimajor axis is the 
tangent plane to the Laplace surface. The transition of the Laplace surface from the ecliptic 
to the equatorial plane occurs approximately at the distance where the torques from the 
solar tides and J2 are equal. Omitting factors of order unity, this distance is often referred 
to as the Laplace radius and is given by Goldreich (1966), 



where Rp and Mp are the planet's radius and mass, and Mq, Op and Cp were previously 
defined. For particles orbiting at large distances from any existing inner satellites, one can 
treat the inner moons' effect as a further contribution to the planetary J2, where (see, e.g.. 




(6) 
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TTN), 

1 

J^Rl ^ J,Rl + o'.rnjM,, (7) 

i=l 

where J2 is the effective J2, and the and rrii are the moons' semimajor axes and masses, 
respectively. Any subsequent references to a planet's J2 in this paper should be understood 
as the effective J2 that includes any inner satellites' contribution to the quadrupole 
potential. 

An important dynamical feature for particle orbits decaying slowly compared to the 
precession timescale (e.g., through P-R drag) is that the orbital inclination to the local 
Laplace plane remains roughly constant (Goldreich 1965). This means that a decaying 
particle orbit starting far from the primary in the planet's orbital plane (which coincides 
here with the local Laplace plane) will have its orbital plane follow the Laplace plane as the 
latter shifts on its way inward toward the planet. 

Following TTN, we neglect any variations in Hq, as well as the precession of rip due to 
the torques on the equatorial bulge from the Sun and other planets. The latter timescale 
is generally much longer than a circumplanetary orbit's precession period, in which case it 
can be safely ignored (Goldreich 1965). 

2.3. The Disturbing Potential 

We now derive the disturbing potential using orbital elements and TTN's notation. 
The obliquity 0q and the variable are defined in Fig. 1, whereas the remaining quantities 
are shown in Fig. 2. 

The eccentricity vector e points toward pericenter and has a magnitude given by the 
orbit's eccentricity. The vector j lies along the orbit normal and has a magnitude chosen 
as Vl — e^, so that j ■ e = and + = 1. We choose to measure the longitude of the 



- 11 - 



z 




n, 



1© 



Fig. 2. — Reference plane (white) is the local Laplace plane. Hp is the planet's spin pole and 



Uq the ecliptic pole. The particle's orbital plane (shaded) is defined by its orbit normal j, 
which can be given in terms of the orbit's inclination i and longitude of the ascending node Q. 
The orientation of the orbit (not shown) within the orbital plane is defined by the so-called 
eccentricity vector, which points toward pericenter e, and is parametrized by the argument 
of pericenter u, measured along the shaded orbital plane, from the line of ascending node. 
We choose to measure the longitude of the ascending node Q from the direction defined by 
Aq X hp. Note that j and e are not unit vectors. 

ascending node Q from the direction along n0 x rip. 



The perturbing potentials (averaged over both the particle and solar orbits) are (see 
Eq. 12 in TTN): 



where \&p and \E'0 represent potentials non-dimensionalized by the factor GMp/a. The 
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quantity e© was defined in Eq. 5 and is given by 

= (9) 

Since our aim is to approach the problem using orbital elements, we employ them to rewrite 
the three scalar products in Eqs. 8. This is perhaps easiest to do by first writing j, rip, 
and e in terms of their Cartesian components. This process yields 

j ■ Hp = (1 — e^)"'^^^(cos0cosi + sin^sinicosw) 

j ■ n0 = (1 — e^)"'^/^[cos (00 — 0) cosz — sin(</)0 — 0) sinzcosfi] (10) 
e ■ n0 = e[sin((;/)0 — (;/))(cosz cos sincu + coscu sinr2) + cos(00 — 0) sini sincj]. 



Substituting into Eqs. 8, we obtain the dimensionless disturbing function R = —{'^p + ^0), 

R=~~l^l TTT, — ^ o\t/o [1 - 3(cos0cosz + sin^sinzcosl])^] + 

8 I 6[1 — e'^yi'^ \ a J 

5e^[sin(00 — 0) (cos i cos sin + cos w sin!]) + cos(00 — 0) sin i sin w]^ — 
(1 — e^)(cos(00 — 0) cosi — sin(00 — 0) sini cosfi)^ — 2e^|', (11) 

where we have written the potential only in terms of e0 in order to explicitly bring out the 
dependence on the semimajor axis a relative to the Laplace radius (note that from Eqs. 
6, 9 and 5, ep/e© = (r^/a)^). 

Implicit in Eq. 11 is a relation between a, and 00, since the semimajor axis 
sets the location of the local Laplace plane and therefore of the z axis (see Fig. 2). 
Changing a therefore alters 0. The transcendental equation that connects these quantities 
is (Tremaine et al. 2009) 



sm . 



tan20= ""^ 5. (12) 



cos 200 + 2yri/a ^ 

Equation 12 has four solutions in a 2tt interval: 0, + vr, and ± 7r/2. The solution 
corresponding to the classical Laplace equilibrium has the property — 00 for a ^ r^. 
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Orbits locally aligned with this surface are stable to small perturbations in their orientation. 
The other two solutions 90° away are always unstable: if the orbit plane is displaced slightly 
from either of these directions, it will precess about the stable solutions; we therefore do 
not discuss them further (see TTN). 

Eq. 11 is equivalent to the disturbing function in Eq. 1.3 of Lidov & Yarskaya (1974), 
which is instead referenced to the ecliptic. The two are simply related by a rotation by 
00 — about the x-axis, where is given by Eq. 12. 

2.4. Dynamics in the Laplace Plane 

While orbits in the classical Laplace plane are always stable against perturbations to 
their orbit normal, they are not always stable to small changes in their eccentricity. We now 
find the radial distance at which the Laplace plane becomes unstable to such eccentricity 
perturbations by considering the limit i — > 0. As the orbital plane approaches the Laplace 
plane, however, Q becomes ill-defined. One therefore typically switches to the variable 
w = Q + u, which smoothly approaches the angle between n© x rip and pericenter (see Fig. 
2). 

Before continuing, we wish to relate the angular variables to ue, since Eq. 2 shows that 
it is specifically the halting oi uje that can create large eccentricities. Careful inspection of 
Fig. 2 shows that orbits in the Laplace plane (the white reference plane) must intersect the 
ecliptic plane (not shown but perpendicular to Hq) along the vector n© x rip. If we make 
the same choice of n© x rip as the arbitrary reference direction in the ecliptic plane, this 
means that orbits in the Laplace plane always satisfy Qe = 0- The angle ue is then just 
the angle from Uq x rip to pericenter, or w. Therefore, for orbits in the Laplace plane, we 
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set z = in Eq. 11 and, writing + a; = w^;, we obtain 



I 3(1 -^e^)3/2 ("^) {^^~^ sin^(</'0-</') sin^ 0JE-{l-e^) cos^ 



(13) 



We can now employ Lagrange's planetary equations to find the orbital elements' time 
evolution. The equation for the eccentricity is (cf. Murray & Dermott 1999, p. 251, noting 
that we have non-dimensionalized R and that over-dots denote time derivatives). 



(14) 



n e dw e Oue 

where we have ignored a term involving the mean longitude at epoch that disappears in our 
orbit-averaged equations. Plugging Eq. 13 into this equation, 

- = ^eoe(l - e^) V2 ,i^2^^ _ ^) 2ue. (15) 
n 8 

Since we are considering orbits in the Laplace plane, j lies along z in Fig. 1, and since is 
refers to the angle between j and the ecliptic axis Hq, tE = (pQ — (p- Hence Eq. 15 matches 
the classical Kozai result (Eq. 2). This is what one would expect as the planet's oblateness 
does not contribute secularly to the eccentricity evolution (e.g., Danby 1962). 

The precession of pericenter is altered by the planet's oblateness, however. The 
appropriate equation is (cf. Murray & Dermott 1999, p. 251), 

— = — = (16) 

n n e oe 

where we have omitted a term oc dR/di that vanishes in the limit i — > 0. Again substituting 

for R from Eq. 13, 

CUE 3ee(l-e2)V2 



n 



|2-5 sin^(00-0) sin^ cue-cos^(0q-0)- _ ^^2^/2 (~) [^-3 cos^ 0] | . 

(17) 

In the limit a ^ r^, </>—)■ 0, and ue (3/2)ep(l — e^)~^ (where, once again from below 
Eq. 11, {rL/aY = Cp/e©)- This matches the rate of precession Wp for the longitude of 



- 15 - 



pericenter relative to the planet's equatorial plane due to oblateness (Danby 1962). In the 
limit a ^ r^, 0o — J- (Eq. 12) and ue — ^ (3/4)eQ(l — e^)-*^/^. Since we have restricted 
ourselves to orbits in the Laplace plane, this does not agree with Eq. 4; rather it provides 
the rate one would obtain solely from solar perturbations after imposing the Laplace-plane 
condition that Qe = (see the paragraph preceding Eq. 13). 

Recall that inducing large eccentricity amplitudes relies on the pericenter precession 
rate approaching zero, i.e., keeping ue constant (cf. Eq. 2). In the above two limiting 
cases, Ue > 0, so if ue is to cross through zero, it must do so at intermediate semimajor 
axes. To find the least stringent condition for the pericenter to lock, one can pick the most 
unstable configuration (i.e., the orientation that generates the largest negative terms). As 
in the Kozai case, this corresponds to siio^ue = 1, i.e., ue = ±90°. Setting cue = in Eq. 
17 with Ue = ±90° yields, to first order in e. 




One cannot analytically find a solution for a since 0, 0q and a are all related through the 
transcendental relation in Eq. 12. But one can see that 4sin^(0Q — 0) will only be large for 
far from 00, and the last term (from J2) is only negative for > 55°. Since is bounded 
to be between zero and the obliquity 00 (see Fig. 2), this suggests that high obliquities 
are required for cj^; to drop below zero. It also means that the roots of Eq. 18 (if they 
exist) should be close to r^, since this is where the Laplace plane transitions and is the only 
situation where (pQ — (f> and r^/a are simultaneously appreciable. 

When Eq. 18 is solved numerically for various 0©, no solution appears for 0© < 68.875°. 
Below this obliquity, orbits are therefore always stable. Beyond this threshold obliquity, 
however, the pericenter can halt for a range in a, thereby generating large-amplitude 
eccentricity oscillations; our value for the critical 00 agrees with that derived differently by 
TTN. 
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2.5. Uranus: A case study 

Uranus is a solar system example with an extreme obliquity {(pQ ~ 98°) beyond the 
threshold value of 68.875°. Hence, circumplanetary particles within a certain semimajor 
axis range will generate large eccentricity values. This unstable range is depicted in Fig. 3, 
which prescribes the minimum non-dimensionalized precession rate Cje/^ as given by Eq. 
17 for low eccentricities and ue = ±90°, plotted vs. semi-major axis (using Eq. 12 to solve 
for (p). 
1 d(U£ 



n dt 




Fig. 3. — For low eccentricities, the minimum UE/n (at uje = ±90° in Eq. 17) as a function 
of semimajor axis. The semimajor axis is in units of the Laplace radius Q4Rp for Uranus 
(Eq. 6). The non-dimensionalized precession rate is expressed as a fraction of the rate for 
tl {CjEln = 3e0/4). In the radial range where {CjE/n)min < 0, tOE/n will cross through 
for certain values of ue- In this radial range, the Laplace plane is unstable to eccentricity 
perturbations. 

Fig. 3 shows that a circular orbit lying in the classical Laplace plane will be unstable 
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in the approximate range 0.93ri < a < l.llri. In the case of Uranus, the effective J2 
including the contribution of the inner satelhtes is approximately 0.019 (Eq. 7), ^ GARp, 
and the unstable range translates to 59.5-Rp < a < 74.9-Rp. 

Fig. 4 displays a numerical integration of a nearly circular orbit (initial eccentricity of 
10"^) started far from the planet in the ecliptic (coincident with the local Laplace plane). 
The particle is then slowly brought inward according to a = age"*/'^, where r = 2.5Myr. 
This interval is much longer than the secular timescale on which the orbit evolves of 
~ O.OlMyr. The functional form for the semimajor-axis decay was chosen to match that for 
P-R drag (Burns et al. 1979); this is simply for consistency with Sec. 6 where we consider 
small particles that are subject to this dissipative force. Uranus' orbit is taken as circular, 
and the whole third body effect of the Sun is included. As our later integrations will include 
radiation pressure, we used the well-established dust integrator dl for all our numerical 
simulations (Hamilton 1993; Hamilton & Krivov 1996; Hamilton 1996; Hamilton & Kriiger 
2008; Tamayo et al. 2011; Jontof-Hutter & Hamilton 2012a,b). 

Since the particle starts in the ecliptic, the inclination relative to the ecliptic tE begins 
at zero. As the orbit approaches the Laplace radius (64 Rp), the inclination follows the 
local Laplace plane toward Uranus' equatorial plane. However, both the eccentricity and 
inclination become unstable immediately upon entering the unstable range at ^ 74.9-Rp. 
One can see this corresponds to the point where ue (bottom panel) stops precessing (at the 
maximally unstable orientation of 270°). The reason the eccentricity does not grow the first 
time oje drops to zero at a ~ 77-Rp is that ue is just above 270°, where according to Eq. 2, 
the eccentricity shrinks. The second time ue < 270°, so e > 0. Once the eccentricity grows, 
e eventually becomes large enough that the second-order eccentricity contribution to the 
last term of Eq. 17 becomes important and causes precession to resume, i.e., the particle 
gets close enough to Uranus at pericenter that J2 re-initiates precession. 
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Fig. 4. — Numerical integration of an initially nearly circular orbit started in the ecliptic at 
120-Rp and slowly brought inward. The top panel plots inclination referenced to the ecliptic, 
so initially is = 0- The middle plot displays the eccentricity history, and the bottom 
plot shows the evolution of uje- The eccentricity and inclination become unstable when 
the semimajor axis reaches ~ 74.9Rp. Note also that this is the point where ue remains 
constant, near the most unstable orientation uje = 270°. 

As mentioned at the beginning of Sec. 2.2, this behavior differs from that of Kozai 
cycles. The oscillations are not regular and the eccentricity and inclination are not coupled. 
In the Kozai case this coupling was due to the conservation of angular momentum along 
the ecliptic axis. The planet's oblateness spoils this symmetry from the Kozai problem 
because it allows the particle to exchange substantial angular momentum with the planet 
at pericenter along that previously conserved direction. Note that when the eccentricity 
becomes large and the particles at pericenter approach the inner moons, the approximation 
used in our calculations and simulations of treating the inner satellites as a contribution 
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to the planet's J2 is no longer an appropriate assumption. However, one would expect 
collisions to remove particles shortly after their orbits cross those of these large satellites. 

We also point out that even though the classical Laplace plane is defined only 
for circular orbits, eccentric Laplace equilibria also exist (TTN). The circular and 
eccentric equilibria bifurcate when the circular solution becomes unstable, and a decaying 
particle-orbit can transfer onto the eccentric-equilibrium track. However, TTN find that 
the eccentric equilibrium becomes unstable almost immediately upon bifurcating from the 
circular solution (see Fig. 6 of TTN). This is especially true for particles starting far from 
the planet and evolving inward through the unstable region (rather than starting close and 
evolving outward). The point at which the classical Laplace plane becomes unstable (found 
above) is therefore a good proxy for when a circular orbit originally in the Laplace plane 
destabilizes. 

3. DYNAMICS IN THREE DIMENSIONS 

By restricting the above discussion to orbits that lie in the local Laplace plane, 
we reduced the dimensionality of the problem to two dimensions. We now address the 
general situation where orbits are inclined to the local Laplace plane. Then the problem 
is inherently three-dimensional, and the enlarged phase space makes it difficult to provide 
detailed general results. Accordingly, we do not pursue a complete analytical theory 
and instead limit ourselves to a qualitative description of orbital behavior based on our 
understanding derived from the above analysis as well as various numerical integrations. 
There are some integrable cases considered by Lidov & Yarskaya (1974); however, in the 
case of interest with finite eccentricity and a ~ r^, these only apply to obliquities of zero 
and ninety degrees, or to polar orbits with the orbital axis pointing along the intersection 
between the planet's orbital and equatorial planes. 
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We can first gain some insight by investigating the equations of motion relative to the 
ecliptic plane. Since J2 perturbations have no secular effect on e, the eccentricity evolution 
depends only on solar perturbations; it is therefore given simply by Eq. 2. In Eqs. 3 and 
4 for di/dt and duE/dt, one would have to add the complicated effect of J2 referenced to 
the ecliptic plane. This is obtained by taking the J2 contribution to R in Eq. 11, i.e. the 
term involving {ri/a)^, setting = 0q, appending a subscript E to all angular variables, 
and applying (cf. Danby 1962, recalling that our disturbing potential is non-dimensional) 

iE 1 f ■ dR dR^ 

— = , A CSC lE^ cot lET, — \ (19) 

n Vl - I OQe dujEJ 

ue VT^dR 1 dR 

— = ^ /. 9 ^g^- 

n e de a/1 - e2 otE 

The resulting equations are complex and difficult to pursue analytically. One can, however, 
gain insight from investigating the effect of a non-zero inclination on the solar perturbations 
that dominate the particle's early evolution far from the planet, before the J2 terms become 
important. 

From Eq. 4, a larger inclination acts to lower duE/dt, bringing the orbit closer to 
instability. This leads to Kozai oscillations for sufficiently large in this limit that ignores 
J2. In this sense, the previous section's situation where a particle begins far from the 
planet in the Laplace plane (where tE = 0) furnishes the best-case scenario for stability; 
duE/dt would have to be substantially decreased by J2 in order for duE/dt to drop to zero. 
One should therefore expect that, if a circular orbit starting in the Laplace plane becomes 
unstable, any orbit initially inclined to the ecliptic will also destabilize. Furthermore, 
inclined orbits should become unstable earlier during their inward evolution than their 
uninclined counterparts would. 

Numerical integrations support this assertion. Figure 5 shows the evolution of nearly 
circular orbits that are slowly evolved inward after beginning far from Uranus at various 
initial inclinations to the ecliptic. One sees that more inclined orbits become unstable 



- 21 - 



earlier in their inward migration. We point out, however, that the orderly progression in 
Fig. 5 results from starting all the integrations with the same initial conditions (other 
than the inclination). In the general 3-D case, all the orbital elements affect the dynamics. 
Altering the initial conditions changes the phase in the elements' evolution at which they 
enter the unstable range in a that we found in the 2-D case; this can change the semimajor 
axis at which the eccentricity grows by ~ 10%. It nevertheless remains true that for a given 
set of initial conditions, increasing the inclination destabilizes the orbit earlier in its inward 
evolution. 




a(Rp) 



Fig. 5. — Orbital eccentricity histories for particles begun far from Uranus (120-Rp) with e = 
10~^ at varying inclinations to the ecliptic. Like Fig. 4, the semimajor axis is brought inward 
according to a = aoe~*^^, with r = 2.5Myr. The figure plots eccentricity vs. semimajor axis, 
where constant offsets have been added to the eccentricities to separate the different plots. 
Higher-inclination orbits are inherently less stable and undergo large- amplitude eccentricity 
oscillations sooner in their inward evolution. 
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Changing the initial e, however, does not have a strong effect on the location at which 
the pericenter halts and the orbit undergoes large-amplitude eccentricity oscillations. This 
can be seen from Eq. 17. A non-zero e enhances the J2 contribution, pushing outward the 
location at which e grows slightly; however, this term's steep dependence on semimajor axis 
of (r^/a)^ allows a small change in a to accommodate a large initial eccentricity. The edge 
of the unstable region therefore shifts by less than a few percent for e < 0.3. 

A more complete investigation of inclined orbits is beyond the scope of this paper. 
We note, however, that the threshold obliquity of 68.875° found by TTN, and derivable 
from Eq. 18, applies to orbits in the Laplace plane, i.e., the most stable configuration. For 
orbits initially inclined to the Laplace plane, the threshold obliquity would be lower. An 
inward-evolving object with an initial inclination close to the threshold value for Kozai 
oscillations (~ 39.2° for low eccentricities) could undergo large-amplitude eccentricity 
oscillations in systems with more modest obliquities. We have verified this, finding that 
for a hypothetical planet with obliquity (/)q = 60°, orbits begun with e = and «s ^ 7° 
undergo eccentricity oscillations in the transition region. The maximum eccentricity 
attained increases with initial inclination, varying from e^ax ~ 0.35 for initial ie = 10° to 
Gmax ~ 0.95 for Ie = 35°. Even at Saturn (00 ^ 27°), initially circular orbits started with 
is = 35° undergo oscillations with a maximum eccentricity of ~ 0.2. 

4. THE EFFECTS OF NON-GRAVITATIONAL FORCES 

The mechanism discussed in this paper occurs when the effects of planetary oblateness 
and solar gravity balance at a ~ r^, (Eq. 6) so as to halt pericenter precession. It is 
therefore appropriate to consider whether additional perturbations might instead keep the 
pericenter moving, thereby stabilizing the orbit. While the previous discussion applies 
to objects of arbitrary size and mass, we now consider radiation forces, which are most 
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important for small dust grains. 

Radiation forces are of particular interest because they provide a natural mechanism 
(P-R drag) for the semimajor axis of dust-grain orbits to decay in toward the planet 
and reach the unstable range in a (Burns et al. 1979). Two further effects that are 
generally important for dust grains are direct solar radiation pressure (Burns et al. 
1979) and electromagnetic forces due to the planetary magnetic field (Hamilton 1993; 
Hamilton & Krivov 1996; Burns et al. 2001). 

At all the solar system's planets except Jupiter, the instability occurs beyond the 
magnetopause, rendering perturbations from the planetary magnetic field irrelevant. 
Furthermore, radiation pressure (discussed below) can remove small particles by pumping 
their eccentricities close to unity. Particles then either crash into the primary or escape the 
system entirely (see Hamilton & Burns 1992). For dust grains starting far from the planet 
(~ 200-Rj,), only particles larger than roughly a few microns in radius survive (see below). 
For particles of this size and larger, even inside the magnetosphere, the planetary magnetic 
field is not important (cf. Fig. 11 in Burns et al. 2001); we therefore ignore it. 



4.1. Radiation Pressure 

Solar radiation pressure, however, can have powerful effects. This perturbation has been 
extensively studied, usually by approximating the planetary orbit as circular and averaging 
over the particle's orbit, which generally changes much faster than the planet's orbital 
period (Burns et al. 1979; Hamilton 1993; Juhasz & Horanyi 1995). Mignard & Henon 
(1984) found an exact solution under these assumptions in a frame rotating with the Sun, 
employing several changes of variables; unfortunately, the inverse transformations to the 
orbital elements that we have utilized are complicated. We therefore choose to instead work 
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in the same inertial system we employed above and to find approximate equations sufficient 
for our needs. 

Upon averaging over tlie particle's orbit, two fundamental timescales remain. The 
ffist is the secular rate at which the orbital elements change (~ nFrad/Fg, where n is the 
particle's mean motion and Frad and Fg are the radiation pressure and planet's gravitational 
forces, respectively); the second is simply the Sun's mean motion about the planet Uq. The 
dynamics are set through their ratio Z = (SnFrad) / {"^nQFg) , where the factor of 3/2 results 
from the exact form of the equations of motion (Burns et al. 1979). Note that since we 
will later be interested in higher-order eccentricity terms, we have removed the changing 
factor of a/1 — from the definition of Burns et al. (1979) so that Z is constant (at a given 
semimajor axis). One can express Z as 



where Qpr is the radiation pressure coefficient averaged over the central star's spectrum, p 
is the particle density and s is the particle radius. Note that smaller particles, with larger 
surface-area-to- volume ratios, are more affected by radiation pressure (i.e., have higher 
Z). However, once particles shrink below the scale of the star's peak emission wavelength, 
they lose the ability to couple to the radiation field and Qpr drops to zero; this occurs at 
~ 0.1 yum for solar- type stars. 

If Z approaches one, radiation pressure pumps a particle's orbital eccentricity to unity, 
most often resulting in collision with the planet. This provides a minimum particle size to 
consider. However, because Z increases with a, this limit would vary with initial location 
from the planet. This is because the importance of radiation pressure relative to the 
dominant planetary gravitational field increases the farther out one orbits in the primary's 
gravitational well. For a more detailed analysis of this threshold and the ultimate fate of 
these grains, see Hamilton & Burns (1992). 




(20) 



- 25 - 



The inclusion of radiation pressure introduces high-frequency variations to all the 
orbital elements at the Sun's orbital rate about the planet, Hq (Burns et al. 1979). Since 
these rates are much faster than the precession rates due to J2 and solar tides, one can 
average over these fast solar oscillations. As shown in the next section, to first order in 
the eccentricity, no secular change in ue occurs. For orbits in the Laplace plane with low 
eccentricities, Eq. 18 therefore remains the averaged condition for the pericenter to halt (at 
Ue = ±90°). Thus, as we argued following Eq. 18, Ue will still only cross through zero at 
approximately the semimajor axis where the Laplace plane transitions from the ecliptic to 
the equatorial plane. However, we will show in Sec. 4.3 that by inducing a slower secular 
change in fi, radiation pressure shifts this transition location. Again using the Uranian 
system as an example. Fig. 6 shows the numerical integration of a Z = 0.1 particle started 
at 80-Rp with a small seed eccentricity and inclination to the ecliptic. The inclination follows 
a modified Laplace plane (cf. Fig. 4) as P-R drag slowly decreases the semimajor axis, and 
Ue only becomes stationary when this transition occurs at a ~ 54i?p. 

4.2. Secular Precession Rates 

Although the equations of motion become more difficult to solve with each added 
perturbation, we can make some analytic progress for small eccentricities and inclinations to 
the ecliptic. The relevant equations of motion, subject to the simplifications mentioned at 
the beginning of this section, are provided by Hamilton (1993). The elements are referenced 
to the ecliptic plane (as there is no ambiguity, we henceforth omit the 'E' subscripts), and 
since we limit ourselves to low inclinations (ignoring terms of order i^), we switch from uj to 
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Fig. 6. — Orbital integration of a Z = 0.1 particle begun at 80Rp from Uranus with e = 10~^ 
and tE ~ 0.06°. Radiation pressure has caused the location of the Laplace plane's transition 
to shift inward from a ~ 75Rp to a ~ 55Rp (cf. Fig. 4). 

the variable w = Q + u. This yields 



n 



—UqZVI — sin(n0t + 5 — w) 
cos[nQt + — to), 



UqZc 



sin(t<7 — Q) sin(n0t + 6 — Q), 



(21) 



where S is the angular location of the Sun at t = relative to the inertial reference direction. 
We approach this system of coupled differential equations through the method of successive 
approximations. Expanding Eqs. 21 in powers of e, we begin by ignoring the terms of order 
e and higher. In this limit, Q is constant, and the solution for the first two equations is 
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given by Burns et al. (1979) in terms of the new variables k = ecosw and h = esinw. 

k = ko — Z cos((5) + Z cos(r20t + 6) 

h = ho- Z sm{5) + Z sin(n0t + 6), (22) 

where in their solution, the time t = has been redefined so that 6 = 0. These solutions 
have a readily-visualized geometric interpretation. The system evolves at a rate Uq along 
the locus of points defined by a circle of radius Z. The center of this circle is offset from 
{ko, ho) away from the Sun's initial position S (see Fig. 7). One can visualize the evolution 
of e and w from such a representation since the orbital eccentricity at a particular point 
in h, k space is given by the distance from the origin, and w by the polar angle. The 
Sun's initial location therefore determines the range of eccentricities explored by setting the 
location of the circle's center. 

We now refine our solution by including terms of order e in an expansion of Eqs. 21 in 
powers of e. Omitting the first equation, which is unchanged, 

117 = Too n^Ze cos{nQt + 6 — w) 

Q = —UqZc sin(cc7 — Q) sin(n0t + 5 — fi), (23) 

where wo is the zeroth-order rate employed in our first solution. Expanding the 
trigonometric functions in the equation for Cl, and using the substitutions k = ecosw and 
h = e sin w, one obtains, 

Cl = —nQZ{hcosQ — A;sinr2)[sin(n0t + 6) cosQ — cos{nQt + 6) sinQ]. (24) 

We now feed our zeroth order solution back into the above equation. In particular, we 
treat Q as constant, and use Eqs. 22 for h and k. Since we are interested in how radiation 
pressure interacts with solar tides and oblateness on long secular timescales, we additionally 
average over a solar cycle from t = to t = 27r /uq. This then yields the simple expression 

<n>^ (25) 
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h = e sin in To Sun 




k = e cos UJ 



(kg - Zcos6, hg - Zsin6) 



Fig. 7. — Geometrical representation of Eqs. 22. The system begins at {ko, Hq) and evolves 
along the perimeter of the circle of radius Z at a constant rate Uq. The orbit's eccentricity e 
and w can be read as the system point's distance from the origin and polar angle, respectively. 
The circle's center relative to (fco, ho) is set by the Sun's initial position, 6, and lies at the 
point {ko — Z cos 6,ho — Z sin 6). 

For small values of Z, this expression is consistent with our previous assumption that Q 
evolves at a rate much slower than Uq and matches numerical integrations well. As Z 
approaches unity, our approximations worsen. 

Applying the same procedure of inserting the zeroth-order solution and averaging over 
a solar cycle in the expression for w in Eq. 23 yields the same value 

< ^ >^ _!^. (26) 

Since w = Q + u, this means that the secular change in w is entirely due to Q. Therefore, 
to our level of approximation, u does not move secularly. This justifies our claim from Sec. 
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4.1 that Eq. 18 still represents an averaged condition for the halting of peri center with 
radiation pressure included. However, as we argue in the next section, radiation pressure 
can change the semimajor axis at which the eccentricity becomes unstable by modifying the 
Laplace surface. 

4.3. The Modified Laplace Surface 

We showed in the previous section that, to first order in e, a; does not move secularly; 
therefore Eq. 18 still holds as a condition for the eccentricity to grow to large values. 
However, the regression of Q in the ecliptic plane induced by radiation pressure (Eq. 25) 
spoils the Laplace equilibrium between solar tides and planetary oblateness given by Eq. 12. 
Radiation pressure creates a modified Laplace surface (on which the torques from all three 
perturbations balance) and shifts the location where the local equilibrium plane transitions 
from the ecliptic to the equatorial plane. 

For prograde particles, and to first order in e and i, the secular regression of the node 
due to solar tides is given hj Q = — (3/4)e0n (e.g., Carruba et al. 2003). Radiation pressure 
thus enhances the nodal regression induced by the Sun. As a result, the semimajor axis 
at which these torques balance those from the planet's oblateness (i.e., the point at which 
the Laplace plane shifts) must move inward, where the effects of the zonal harmonics are 
stronger. 

We do not calculate the detailed warp of the Laplace plane in this paper — for references 
on the process in the classical case of the competition between solar gravity and planetary 
oblateness, see Ward (1981) and Dobrovolskis (1993); Allan & Cook (1964) treat the 
general case of an arbitrary number of non-interacting perturbers in independent planes; 
for the different case of radiation pressure offsetting planetary oblateness, see (Hamilton 
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1996, cf. Fig. 18 of Burns et al. 2001). In our case involving all three perturbations, we 
limit ourselves to estimating the transition location of the Laplace plane where we expect 
the orbit to become unstable and execute large-amplitude eccentricity oscillations. More 
concretely, this will approximately correspond to the distance at which the nodal precession 
rates from the solar tides and radiation pressure balance that due to planetary oblateness: 
^Sun + ^Rad = ^J2- One cau understand this as an approximate condition that the torques 
from forces directed out of the orbital plane cancel (see Eq. 38 in Burns 1976). The nodal 
rate due to J2, now referenced to the planet's equatorial plane and again to first order in 
the inclination, is Qj2 = — (3/2)epn (e.g., Murray & Dermott 1999, p. 270). The nodal rate 
due to solar gravity is flsun = — (3/4)e0n (e.g., Carruba et al. 2003). 

Substituting for e0 and ep from Eqs. 5 and 9, the condition for the balance of precession 
rates is 



where is the approximate semimajor axis at which the Laplace plane transitions. For 
Z = 0, the equation can be solved analytically, yielding qt = 2^/^ri. At Uranus, this 
corresponds to 74.7i?p, which can be seen from Fig. 4 to be approximately the location 
where the inclination is intermediate between and Uranus' obliquity of 98°. It furthermore 
quite accurately matches the location at which the eccentricity grows rapidly, and is 
therefore a slightly more accurate estimator than the dimensionally obtained ri. 

One can roughly estimate the particle size-range in which radiation pressure is important 
at the Laplace plane transition by setting Clrad ~ ^Sun, or nQZ{r-iY ~ ^Q'^ijL) = nQ/n{ri), 
where is the Laplace radius from Eq. 6. This will generally yield a small value of Z 
since Z = 1 would correspond to a Laplace radius equal roughly to the Hill radius. In the 
Uranian example previously discussed, this corresponds to Z ~ 0.05, or for particles with a 
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density of Ig/cm^, to s ~ 70/im. Figure 8 shows numerical integrations of particles around 
Uranus with radii s = 20/im-80/im, with the corresponding a-r (numerically obtained) 
marked vertical line. 
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Fig. 8. — Orbital integrations of particles with various radii s orbiting Uranus. The four 
panels, from top to bottom, correspond to values of Z (at a = 75Rp) of 0.04, 0.07, 0.11, and 
0.17. Particles were started at a = 90Rp with a seed eccentricity and inclination of e = 10"^ 
and i = 0.06°, respectively. The vertical solid lines denote the transition locations of the 
Laplace plane for each size computed from Eq. 27. The dashed lines denote the transition 
location in the absence of radiation pressure. This predicted position matches the location 
where the eccentricity destabilizes to within ^ 10%. 



The modified location of the Laplace-plane's transition from Eq. 27 matches the 
onset of instability to within ^ 10%. As the particle size decreases and Z increases, our 
approximations worsen, and one can see in the fourth panel that the behavior is beginning 
to change; the eccentricity first decreases, and later rises in two steps. Our results should 
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therefore be applied with caution beyond Z > 0.2. We find that for large values of Z, some 
particles retain low orbital eccentricities as they traverse the unstable region. We note, 
however, that this range in Z represents a narrow size range since Z oc . In this example, 
the range Z = 0.2 — 1 only corresponds to s ~ 20/im-4yU,m (particles with Z > 1 need not be 
considered as they would have been immediately removed). If interested in these smallest 
particles, one must carry out suites of numerical integrations over a wide range of initial 
conditions to capture the full dynamics. 



4.4. Retrograde Orbits 

We now briefiy consider retrograde orbits, which interestingly can exhibit qualitatively 
different behavior. For retrograde orbits, Qsun and ^lj2, which both contain a cosz 
dependence, switch sign. One can obtain Clrad by rederiving the results of Sec. 4.2 starting 
from the equations given by Hamilton (1993) with i ^ 180° instead of Eqs. 21; alternatively, 
one can employ a symmetry argument similar to ones presented in Hamilton (1994). 

One can change a retrograde orbit into a prograde orbit by rotating the coordinate 
system by 180° around the x axis, so that z — )■ —z. One can then immediately write 
down the solution found above for prograde orbits, except in this coordinate system the 
Sun now moves retrograde, so one must make the transformation Uq — )■ —Uq. This yields 
^rad ~ +'^0^^/2j where the superscript minus sign denotes that these are elements in the 
—z coordinate system. The final step is to relate to Q~^, the longitude of the ascending 
node in the original coordinate system. Since, by the right-hand-rule, the directions in 
which angles increase in the +z and —z coordinate systems are opposite in direction, 

= —Q^. That actually is not quite right, since upon flipping the conventional "up" 
direction, the ascending and descending nodes switched places, so fi"*" = 180 — . This 
yields = —Q~ = — ?t,qZ^/2; therefore, while the rates due to solar tides and planetary 
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oblateness flip sign for retrograde orbits, the rate due to radiation pressure does not. This 
is because, while planetary oblateness and solar tides (after averaging over a solar orbit and 
smearing the Sun's mass into a ring) are symmetric under z — )■ — z, radiation pressure is 
not, due to the Sun's motion changing sense. 

The condition for the three torques to balance therefore becomes ll^sunl ~ \ ^rad\ = \^J2\- 
In this case radiation pressure detracts from the solar rate, so the transition location will 
move outward, where weaker oblateness perturbations are sufficient to offset the reduced 
combination. There is the further possibility that \flrad\ overwhelms l^sunl, in which case 
the balance condition cannot be satisfied. Since Ifi^unl oc a'^^^ while \(lrad\ a, there will 
always exist an a at which the solar rate overtakes the rate due to radiation pressure; 
however, if that a lies beyond the particle's initial semimajor axis (which is constrained 
to be smaller than the Hill radius), radiation pressure will always dominate. In this case 
there is no Laplace equilibrium and the inclination does not transition to the equatorial 
plane. The instability is thereby avoided. The threshold Z where this occurs is given by 
the condition \Qrad\ = 1^5™!- We considered the balance of these two rates at the Laplace 
radius in the prograde case; the result evaluated at yields the threshold value of Z, 
Zt ~ [3nQ/(2n)]^/^. Using Eq. 20 to solve for the threshold particle size St, 

The threshold size in the Uranian case with oq = 140-Rp is St ~ 46/im. Fig. 9 shows the 
range of behaviors discussed in the previous paragraph for the Uranian case with the same 
particle sizes as in Fig. 8. The direct integrations match our analytic predictions well for 
our low chosen values of Z. The irregular behavior in the 50/im case is presumably the 
result of its proximity to the threshold size from Eq. 28, but we do not investigate this 
further in this paper. 
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Fig. 9. — Orbital integrations of retrograde particles with various radii s orbiting Uranus. 
The four panels, from top to bottom, correspond to values of Z (at a = 75Rp) of 0.04, 0.07, 
0.11, and 0.17. Particles were started at a = 140i?p with a seed eccentricity and inclination 
of e = 10~^ and i = 179.91°, respectively. The vertical solid lines denote the transition 
locations of the Laplace plane for each size computed from the appropriate condition for 
retrograde orbits discussed in the text. The dashed line denotes the transition location in 
the absence of radiation pressure. For the bottom two panels, the transition locations are 
at a = 751-Rp and a = 3803i?p, the latter of which is beyond the Hill sphere. In these two 
cases, the Laplace plane does not transition to the equatorial plane and the eccentricities 
remain stable. 

5. CONCLUSION 



We have shown that the unstable range in semimajor axis around planets with high 
obliquities found by Tremaine et al. (2009) can be understood as a modification of Kozai 
oscillations. Furthermore, we extended their work (which focused on orbits lying in the 
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local Laplace plane) and provided equations valid for arbitrary inclination. Although it 
is difficult to give precise general results, we showed that orbits with inclinations off the 
Laplace plane are less stable. We therefore argued that the threshold obliquity of 68.875° 
found by Tremaine et al. (2009) is an upper limit — inclined orbits can become unstable 
around planets with lower obliquities. 

We then investigated the instability as it applies to dust grains. Dust grains are subject 
to Poynting- Robert son drag, which provides a natural mechanism to sweep the semimajor 
axis inward toward the unstable region. However, one must also consider the additional 
effects of radiation pressure on dust-particle orbits. We found that radiation pressure 
modifies the classical Laplace surface, and that this shifts the unstable range of semimajor 
axis. For prograde particles, this chaotic region is shifted inward, while for retrograde 
particles it is shifted outward, and can even disappear for small particles. We estimated 
the threshold grain size at which orbital eccentricities remain stable for retrograde particles 
in Eq. 28. For the smallest particles with Z > 0.2 (cf. Eq. 20), or particles with large 
initial inclinations or eccentricities, our analytical approximations break down. We found 
in simulations that in such cases, for a minority of initial conditions, even prograde orbits 
can remain stable. Suites of numerical simulations spanning the range of initial conditions 
are therefore required to fully characterize a population of dust evolving in toward a 
high-obliquity planet. 

This work can be applied both in the solar system and beyond. Bottke et al. (2010) 
have proposed that, at each of the giant planets, a vast supply of dust generated by 
the irregular satellites once existed. At least in the case of Saturn, this supply persists 
today (Verbiscer et al. 2009). Many irregular satellites have inclinations close to the 
low-eccentricity threshold for Kozai oscillations, i ~ 39.2° or 150.8°. These orbits are very 
unstable, and dust originating from such objects might undergo large-amplitude eccentricity 
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oscillations even around planets with modest obliquities. At Uranus, all but the smallest 
particles will do so, and this might explain the color dichotomies common to the outer four 
regular satellites observed by Buratti & Mosher (1991). Tamayo et al. (2012) have started 
toward such an explanation, which we will pursue elsewhere. More generally, this instability 
could be applied to myriad classes of circumstellar binary objects, such as binary KBOs and 
asteroids. Finally, having incorporated radiation forces, one could consider debris disks in 
systems with an interior planet (providing an effective J2) and a highly-inclined companion. 
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